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SUMMARY 

Two computer programs have been developed for subsonic 
inviscid analysis and design. The first solves arbitrary mixed 
analysis-design problems for multi-element airfoils in two- 
dimensional flow. The second calculates the pressure distribu- 
tion for arbitrary lifting or nonlifting three-dimensional con- 
figurations . 

In each program, inviscid flow is modelled by using 
distributed source-doublet (or source-vortex) singularities on 
configuration surface panels. The prescribed normal velocity 
distribution is satisfied indirectly by applying an internal 
perturbation potential boundary condition to the center of each 
panel. The method is numerically stable, and the prediction 
accuracy is competitive with more complex curved panel formula- 
tions . 


In the two-dimensional mixed analysis-design program, the 
geometry of each airfoil element is divided into analysis regions 
of prescribed geometry and design regions of prescribed pressure. 
Iteration is used to predict both the shape and length of the 
geometry in design regions. Satisfactory convergence is usually 
obtained in only five iteration cycles because the inverse step 
uses the exact derivative matrix of velocity with respect to 
geometry perturbation. Example solutions illustrate applications 
of the mixed analysis-design feature to high lift airfoils and 
viscous-inviscid interactions. Also demonstrated is an improved 
Kutta condition model for airfoils with finite thickness trailing 
edges . 

The three-dimensional program uses constant source and quad- 
ratic doublet distributions on arbitrary flat quadrilateral panels 
to induce the flow field. The method of least squares is applied 
to minimize doublet discontinuities at the panel edges, thereby 
reducing the sensitivity of the numerical solution to the panel- 
ling selected by the user. 

A first-order expansion approach is defined for generating 
the matrix of velocity perturbations with respect to surface 
geometry perturbations. In the three-dimensional program, such 
an approach would allow the user to predict the effect of 
arbitrary configuration changes by the inexpensive process of 
multiplying a precalculated matrix by a new right-hand-side. A 
method of using the matrix for three-dimensional wing design is 
also described. 

Mathematical formulations and representative solutions are 
presented for the computer programs. 
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INTRODUCTION 


The surface panel method philosophy for solving arbitrary 
subsonic potential flow problems involves the mating of classical 
potential theory with contemporary numerical techniques. 

Classical theory is used to reduce an arbitrary flow problem to 
a surface integral equation relating boundary conditions to an 
unknown singularity distribution. The contemporary numerical 
techniques are then used to calculate an approximate solution to 
the integral equation. This involves representing flow boundaries 
by surface panels on which potential flow singularities are 
distributed . 

Whereas there is no limit to the number of different singu- 
larity distributions that can induce a given flow field, the 
type of singularity plays an important role in determining the 
success of a numerical solution method. The advantages of using 
the combined source-doublet distribution corresponding to the 
classical third identity of Green were described and demonstrated 
in reference 1. The combined distribution is not subject to the 
prediction errors associated with the application of source 
methods to thin, highly loaded wings or nacelles. In iterative 
design methods, the combined source-doublet distribution can be 
used to generate converged solution geometries corresponding to 
prescribed pressure distributions. This contrasts sharply with 
doublet (or vortex) methods, which are numerically unstable for 
leading edge design. If indirect internal perturbation potential 
boundary conditions are applied to constant source-quadratic 
doublet distributions on flat panels, the resulting prediction 
accuracy is typically competitive with more complex curved panel 
formulations . 

This report describes recent developments at McDonnell 
Aircraft Company (MCAIR) of two and three-dimensional subsonic 
analysis and design methods using combined source-doublet panels 
to represent the inviscid flow. Formulations are presented for 
two potential flow computer programs. The first solves arbitrary 
mixed analysis-design problems in two dimensions for multi-element 
airfoils. The second solves three-dimensional analysis problems 
for wing-fuselage configurations. 

The mixed analysis-design multi-element airfoil program was 
created by extending an earlier potential flow program that would 
solve all-analysis or all-design problems (reference 1) . In 
regions selected by the user for design, the geometry correspond- 
ing to a prescribed pressure distribution is calculated by itera- 
tion. Both the shape and length of the region will converge on 
the geometry that most nearly corresponds to the prescribed 
pressures. For airfoils with thick trailing edges, a method is 
presented in which the predicted lift is improved by panelling 
the base. 
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The three-dimensional panel program is based on the method 
described in reference 1. An improvement that has been incor- 
porated applies a spline fit technique similar to one developed 
by Boeing (reference 18) in order to reduce singularity dis- 
continuities at panel edges . 

An approach is discussed for generating velocity deriva- 
tives with respect to arbitrary geometry perturbations of a wing- 
fuselage configuration. Application of this approach to the 
three-dimensional panel program would enable the prediction of 
the effect of arbitrary configuration changes by simply multiply- 
ing a pre-calculated matrix by a new right-hand-side. Also 
discussed is a procedure for using the velocity derivatives to 
solve inviscid wing design problems. 

This report presents the mathematical formulation and 
example solutions for the computer programs. 

AIRFOIL MIXED ANALYSIS-DESIGN METHOD 

A variety of two-dimensional incompressible potential flow 
problems can be solved by the Multi-Element Airfoil Inviscid 
Analysis and Design Program (Version 1) , hereinafter designated 
"Program MAAD" . Figure 1 illustrates the type of mixed analysis- 
design problem that can be solved. The geometry of one or 
more airfoil elements is divided into analysis regions of pre- 
scribed geometry and design regions of prescribed pressure dis- 
tribution. It is not necessary to have regions of both type on 
an airfoil; therefore, analysis-only or design-only problems 
are permitted. For each analysis region, a normal velocity 
distribution is prescribed, and the tangential distribution is 
calculated. For each design region, desired normal and tan- 
gential velocity distributions in terms of surface distance are 
prescribed simultaneously, and the corresponding boundary 
geometry is calculated. In the typical design problem, the ob- 
jective is to calculate a flow streamline of given pressure 
distribution. This corresponds to the prescription of zero 
normal velocity, with the prescribed tangential velocity distri- 
bution being generated from the pressure distribution by the 
application of Bernoulli's equation. 

Program MAAD is an extension of an earlier two-dimensional 
multi-element airfoil program described in references 1 and 2. 

The analytical portion of the extension was accomplished as part 
of the MCAIR Independent Research and Development Program. 

Coding, check-out, and evaluation were performed as part of this 
contract. The significant new options are: 

(1) Arbitrary mixtures of analysis and design regions on 
each airfoil element, 

(2) Non- zero normal velocity boundary conditions. 
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Figure 1 . Illustration of a Mixed Analysis - Design Problem 

(3) Sharp corners at arbitrary points on an airfoil (such 
as a flap hinge line) , 

(4) Arbitrary free stream velocity (including zero), and 

(5) Simulation of two-dimensional wind tunnel wall in- 
terference . 

Flow fields are modelled in Program MAAD by a combined 
source-vortex surface panelling method. All surface panel 
methods are exact in the sense that the difference between the 
calculated and analytical incompressible potential flow solutions 
can be made arbitrarily small at the expense of increased com- 
puting time. However, significant differences between the 
various methods exist with respect to accuracy for a given number 
of panels, numerical sensitivity to certain boundary shapes or 
user generated panelling anomalies, computing efficiency, 
applicability to design problems/ and formulation simplicity. 

For two-dimensional flow, nearly all the possible approaches 
depicted in figure 2 were evaluated at MCAIR prior to this con- 
tract. For a variety of geometries, it was consistently 
demonstrated that the approach selected is competitive with or 
superior to each of the others . 

The approach is to solve for the combined source-vortex 
distribution associated with the classical third identity of 
Green (reference 3) . A detailed description of the theory and 
discussion of the advantages of this surface singularity panel 
method are provided in reference 1 and will be summarized only 
briefly here. The reason for selecting combined source-vortex 
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Singularity Type Solid Body Boundary Conditions 

• Source • Zero Normal Velocity 

• Vortex (Doublet) • Constant Stream Function 



Figure 2. Menu of Surface Singularity Panel Methods 

Incompressible, Potential Flow 


distributions is to eliminate the numerical instabilities asso- 
ciated with: (1) the application of source-only panel methods to 

the analysis of thin, highly loaded geometries and (2) the design 
of leading edge regions by vortex-only methods. The use of 
indirect internal perturbation potential boundary conditions, 
flat surface panels, and the local singularity-velocity relation- 
ships of Green's identity lead to an efficient, straightforward 
panel method. Furthermore, the subsequent prediction accuracy 
is typically competitive with more complex higher order (curved 
panel) formulations (references 4 and 5) . It is noteworthy that 
the first use of indirect internal potential boundary conditions 
was by Morino (reference 6) . 

In Program MAAD , the vortex distribution is linear on each 
panel. The program has an option for either piecewise constant 
or linear source distributions. In almost all cases, the type of 
source distribution selected has no appreciable effect on the 
solution . 

Mixed analysis-design and design-only problems are solved by 
iteration. The user must provide a starting geometry for each 
design region in order to initialize the calculations. As illus- 
trated in figure 1, continuity with an adjacent analysis region 
is not required at the aft end of a design region starting geometry. 
Before any aerodynamic calculations are performed, the program 
automatically satisfies continuity requirements by generating the 
least possible stretching and shape perturbation of the panelled 
design region starting geometry. Then the combined source-vortex 
representation is used to calculate an analysis-only solution. 

Next, a first-order inverse method is applied to predict the 
design region geometry modification that will minimize the dif- 
ference between the prescribed and the calculated tangential 
velocity distributions. The inverse method satisfies normal 
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velocity boundary conditions and geometric continuity as first- 
order constraints. Finally, the predicted geometry modification 
is applied in order to provide the starting geometry for the next 
iteration cycle. 

This section presents the detailed mathematical formulation 
for (1) the analysis-only method and (2) the first-order inverse 
method for design regions. It is assumed that the reader is 
familiar with surface singularity theory. Reference 1 provides 
adequate background. The approach used in the first-order in- 
verse method was created in 1973 (reference 7) . Numerical 
examples demonstrating the capabilities of Program MAAD follow 
the formulation presentation. 

Analysis-Only Solution Method 

Given one or more airfoil elements with prescribed normal 
velocity (Neumann) boundary conditions, the objective is to 
calculate the tangential surface velocity distribution in in- 
compressible potential flow. The solution process consists of 
the following five steps: 

(1) Establish a panelled surface geometry. 

(2) Calculate the source distribution based on the local 
surface slope with respect to the free stream 
direction. 

(3) Calculate source and vortex potential influence coef- 
ficients . 

(4) Set up and solve the system of linear equations re- 
lating boundary conditions to unknown vortex strengths. 

(5) Calculate the tangential velocity distribution from 
the solution vortex strengths. 

These steps are described in detail below. 

The geometry panelling nomenclature is defined by figure 3. 
Each airfoil element is identified by index q (1 <_ q <_ q>p) . The 
geometry of an element is defined by n surface points input in 
clockwise order. The trailing edge is defined by points 1 and 
n. Normally, these two points coincide to form a closed panelled 
geometry. Although it is not recommended, a gap between points 
1 and n is permitted. Completely enclosed flows are allowed. 

For example, consider the wind tunnel test section of figure 4. 
The rectangular enclosure is treated as an airfoil element. 

The only difference is that the input points of an enclosing 
boundary are ordered counterclockwise. The velocity components 
at the midpoint of any panel i on any element q are designated 
VTi and VNf. It is noted that positive is directed toward 

the real fluid side of a panel. On the other side of a panel, 
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Figure 3. Multi-Element Airfoil Panelling Nomenclature 



Figure 4. Simulation of Wind Tunnel Wall Interference 
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There is one boundary condition control point at the mid- 
point of each panel. The control point is slightly withdrawn to 
the imaginary flow (internal) side of the panel in order to 
avoid a discontinuity in potential at the exact panel surface. 

At each sharp corner, an additional internal control point is 
inserted (figure 5) to control the discontinuity in vortex den- 
sity that is theoretically expected at a surface slope discon- 
tinuity. Here a sharp corner refers to the intersection of two 
panels at a point where the user designates that the actual 
airfoil surface slope is discontinuous. At a trailing edge, two 
Kutta condition control points are inserted in the actual flow 
to control the circulation (figure 5). It is not necessary for 
a trailing edge to be a sharp corner. Therefore, a smooth geome- 
try such as an ellipse can be modelled. Definitions for the 
terms depicted in figure 5 follow: 


-y -y 



As 


C 



(As a + As b ) 


( 2 ) 


as te = / (as a + As b } (3) 

where As A and Asb are the lengths of panels A and B, and where 
ki and k 2 are input values much smaller than 1. The recommended 
values are K i = 0.01 and K 2 = 0.02. The value for e (figure 5) 
has been arbitrarily selected to be 1/2 radian in Program MAAD . 



Figure 5. Control Points at a Sharp Corner and/or Trailing Edge 
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Figure 6 illustrates the linear source and vortex distribu- 
tions on the panel designated B, adjacent to panels A and C. 

Here the £-axis measures surface distance, so that panels A and 
C can be assumed to have been rotated to coincide with the plane 
of panel B. The source function a(£), which passes through °b» 
is a close linear curve fit to the source density at the three 
panel midpoints (° A , ° c ) . The linear vortex function y(£) 

is determined by the vortex density at the endpoints of panel B 
( Y AB' y BC ) • 


a(C) = 


B 


2 ? 


As + As 

( AS + AS J (a B " °A ) 
A B 


AS A + As b 

+ l As + As * ( °C " °B ) 

rJ U 


AS a + 2As b + AS C 


(4) 




— ( Y tY ) 
2 ky AB r BC' 
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- Y 


AB 
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) Z 
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Figure 6. Singularity Density Distributions on a Panel 
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Special handling is required at sharp corners. Suppose 
that the intersection of panels A and B in figure 6 is designated 
a sharp corner by the user. Then theoretically neither the 
source nor vortex densities will be continuous at the inter- 
section. To accommodate this, a(£) on panel B is selected as 
the linear function passing through ag and ct c , instead of the 
function defined by equation (4) . For the vortex density, yab 
is assumed to be double-valued, i.e., yab an< ^ YAB respectively 
refer to the vortex density immediately to the left and right 
of the sharp corner. Analogous treatment is applied if the 
intersection of panels B and C is designated a sharp corner. 

In accordance with Green's third identity, the source 
strength at the midpoint of each panel i on element q is imme- 
diately determined from the following equation: 


CT 


i 


V >T + V sin ( 6 . - a) 
N . «> 1 


(6) 


where Via. is 

N i P 

equations (4) 
multi-element 


the prescribed normal velocity. By combining 

and (6), the source distribution over the entire 
airfoil system is uniquely determined. 


The potential induced at each control point by the source 
and vortex distributions on each panel is calculated as follows. 
First, a (C , n ) -coordinate system is generated for the influ- 
encing panel as illustrated in figure 7. Then the coordinates 
of the control point of interest are converted from the (x,y)- 
system to the (C , n) -system. Using subscript P to designate the 
control point location. 
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The following functions provide the potential at point (£ p , n p ) 
induced by unit strength and unit linear strength source and 
vortex distributions on the influencing panel. 

Definitions : 


0 (n p i. 0) 

e - 

( (rip < 0) 


(9) 
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Induced Potential: 


Unit source cr{£) = 1 
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By combining equations (4) and (5) with equations (19) -(22), it 
is possible to determine the potential <j> at (?p,np) induced by 
the simultaneous action of all singularities on the influencing 
panel. To accomplish this, it is convenient to cast equation 
(4) in the following form: 


a (?) = 


B 


(< Va + °2°B + °3°C 1 <4i 
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(23) 


Here panel B of figure 6 is considered to be the influencing 
panel. It follows that 
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(24) 


The coefficients of a^, a B , oq, yab , and ybC are designated 
potential influence coefficients. With the exception of yab and 
Ybc ’ all the terms in equation (24) have been previously defined 
as functions only of V K , a, and the input coordinates (x^,y^). 
Therefore, equation (24) can be cast in the following form: 
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(25) 


- d l + 


d 2 Y AB 


d 3 Y BC 


where only y, E and y__ are unknown. 

Aid r>C 

The internal boundary condition to be satisfied is that 
the imaginary internal flow of each airfoil element be of uniform 
velocity V^. This is accomplished by requiring that the internal 
perturbation velocity line integral between any two adjacent 
internal control points i cp and i C p+l be zero, where the control 
points are ordered in the direction of increasing panel index i. 
The line integral is defined as 


I(i CP> s 


i +1 

C 

CP 


(V - v m ) 


ds 


(26) 


where V is the total internal velocity and ds follows an inter- 
nal path, i.e., a path in the imaginary flow. With the aid of 
equation (25), the contribution of all singularities on influ- 
encing panel B to line integral I(i cp ) is expressed as 


I(l CP ) B + (d l + d 2 Y AB + d 3 Y BC ) i Cp +l 


~ (d l + d 2 Y AB + d 3 Y BC ) i 


CP 


+ A 2 (y AB + Y BC ) AS B 


(27) 


where A is equal to 0, +1, or -1, depending upon whether the 
internal path between control points i^p and igp+1 crosses the 
positive 5-axis of panel B. The A-term essentially cancels out 
the discontinuity in vortex-induced potential across a slit. 

If there are an even number of crossings, A = 0. If there are 
an odd number of crossings, A = + 1 if n at point i(^p is < 0. 
For simplicity, the internal path selected is the one that fol- 
lows the two adjacent panels on which i^p and i^p+1 lie. It is 
noted that if there are a total of n C p internal control points 
for element q, the line integral between control points n cp and 
1 is not to be set to zero explicitly. This follows from the 
fact that any closed internal velocity line integral will be 
zero; therefore, there can be only n^p-l linearly independent 
internal boundary conditions per airfoil element. 

It is now possible to establish the system of linear 
equations relating boundary condition satisfaction to the un- 
known vortex densities. The equations to be solved are ordered 
as follows: 
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Now consider arbitrary panel j to be the influencing panel B. 
The unknown vortex density Yab identified as y-; . In the 
event that the intersection between panels A and B is a sharp 
corner , 


y ab Y j 


Y AB " Y j + AY 


( 28 ) 


where Ay is an additional unknown. There is one unknown vortex 
discontinuity Ay per sharp corner. 


The unknown vortex densities are set up as an array u • , 
ordered as follows: 


u 


1 


u 


2 


Y 


1 


▼ 


(n-1) 


Ay- 

I' 

Ay 


(n cp +l-n) 


I 


Element 1 


Element q 


* 


t 


I 


Element q m 
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The total number of equations and unknowns are the same. That 
number is equal to the sum of the number of panels and sharp 
corners in the entire multi-element airfoil system. It remains 
to define the circulation control equation required for each 
element q. 

Any one of five available circulation options can be 
selected for each airfoil element. The basic Kutta condition 
will be described fully; the other four options will be summar- 
ized. The Kutta condition is satisfied by requiring that the 
average velocity between the two control points in figure 5 be 
zero. This is accomplished in the same manner described above 
for two adjacent internal contol points, except that the line 
integral between the Kutta condition control points includes the 
free stream contribution to velocity. The free stream contri- 
bution to the line integral between the upper and lower points 
is V cocos a (xupr - ^lWr) + Vasina (y upR - y LWR ) • If the angle 2e 
between the two control points were allowed to approach zero, 
this Kutta condition would exactly approach the condition of zero 
velocity normal to the trailing edge bisector at the distance 
AsipE downstream. The finite value e = 1/2 provides a close 
approximation to zero normal velocity, without requiring the 
calculation of velocity influence coefficients. 

Following are summaries of the other four circulation op- 
tions ; 

(1) The extrapolated bisector Kutta condition is identical 
to the basic condition described above, except that the 
angle 9 q of figure 5 is replaced by a trailing edge 
angle determined by extrapolation from two upper 

and two lower surface upstream panels. The extrapo- 
lated Kutta condition is recommended over the basic 
Kutta condition. 

(2) A numerical value of circulation can be prescribed, 
which will be set equal to the net integrated vortic- 
ity around the panelled element. 

(3) Circulation normalized by the panelled perimeter of 
the element can be prescribed, which is equivalent 
to prescribing the average vortex density. The only 
practical difference between this option and option 
(2) is for design problems, where the perimeter may 
not be known beforehand. 

(4) If the present analysis solution is part of an iter- 
ation cycle in a design solution and if there are any 
design regions on the element in question, then the 
prescribed midpoint tangential velocity components 

Vm . can be used to fix the circulation. The corres- 
Tl P . 

pondrng crrculation control equatron is 
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E (V, 


T, 


As^) — £ ( Vip _ 


As i ) 


(29) 


where the summation includes design regions only. 

In the above options (2) , (3) , and (4) , the Kutta condition 

control points are ignored. 

If the entire flow is bounded by one element of counter- 
clockwise ordering, the vorticity of that element is fixed by the 
condition that the total vorticity of the q^ elements must be 
zero. This is a consequence of the fact that the imaginary flow 
of the bounding (counterclockwise) element must have zero circu- 
lation in order to satisfy Green's third identity. 

The system of linear equations to be solved is set up in 
the following form: 


2 C . . u . = RHS . (30) 

j 13 3 i 

where Uj is the array of unknown vortex strengths and RHS. is the 
known right-hand-side of the ith equation. A standard matrix 
inverse method is used to solve the system. The matrix inverse 
and potential influence coefficients are stored for future use if 
the next analysis case involves changes in only angle of attack, 
circulation, and normal velocity distribution. It is noted that 
although any such change will alter the right-hand-side RHSi, 
neither the potential influence coefficients nor matrix Cj_j will 
be affected. 

Corresponding to Green's third identity, the solution 
tangential velocity at the midpoint of panel i of element q is 

V T± = V oo cos(0 i - «> + \ + Y (i+1) ] (3D 


If point i is a sharp corner, then y ■ is to be replaced by 
y. + Ay, corresponding to equation (28). If (i+1) = n, then 
y is to be replaced by y^. 

Although prediction accuracy is usually better at panel 
midpoints, the velocity at panel endpoints is also calculated. 

In design problems, both midpoint and endpoint tangential velocity 
components are prescribed to provide numerical stability in the 
inverse step. The tangential component of velocity at point i 
is calculated by the formula 



V cos(0„ - a) + y ■ 

co hj . l 

1 


(32) 
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( 33 ) 


where subscript E designates a panel endpoint, i.e., point i. 
The angle 0 E . is a weighted average of the two adjacent panel 
angles. 


E i 


AS . 6 . . + As . . 0 . 

1 1-1 1-1 1 

AS . + AS . -| 

1 1-1 


If e i-i and 0 i differ by more than ir , it is necessary to adjust 
one of the angles by adding the proper multiple of 2 it before ap- 
plying equation (33) . If i = 1 or n in equations (32) and (33) , 
replace subscript (i - 1) by (n - 1) and i by 1 . 

The only remaining step in the analysis-only solution for- 
mulation is the calculation of pressures, forces, and moments. 
Under the presumption of steady state flow, Bernoulli's law is 
used to calculate the midpoint pressure coefficient for each 
panel i. 




+ V 


N 


V 


ref 


2 


(34) 


where V r ef is usually , unless V m is zero. Force and moment 
coefficients are calculated by the approximations that the 
panelled geometry is the actual surface and that the pressure 
coefficient on panel i is uniform with value Cp^. 

First-Order Inverse 

In design problem solutions, each iteration cycle is divided 
into an analysis-only step and a first-order inverse correction 
to the geometry. The former was described above. The formu- 
lation for the latter is presented here. The objective is to 
determine the change in surface geometry most nearly correspond- 
ing to the prescribed tangential velocity distribution in design 
regions. The approach is to establish to first-order the rela- 
tionship between velocity and geometry perturbations. This re- 
lationship is then inverted to determine the geometry pertur- 
bation corresponding to the desired velocity perturbation. 

Before proceeding with the solution formulation, it is 
necessary to precisely define analysis and design regions. The 
indexing nomenclature for the division of an airfoil element 
into regions is presented in figure 8. Region p is bounded by 
i A (p) i <_ i B (p), where i is a panel endpoint index defined 
earlier in figure 3. At this stage, it is assumed that all 
continuity requirements are satisfied. For p ={= p T , the con- 
tinuity requirement is that points i B (p) and i A (p+l) coincide. 

For p = prp, the requirement is that the components of the trail- 
ing edge opening (x n -xq, y n - Yl) equal prescribed values. Table 
I summarizes the differences between analysis and design regions. 
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Figure 8. Division of an Airfoil Element into Regions 


Region Type 

Boundary Conditions 

Geometry (S ha 

pe and Length) 

Normal 

Velocity 

Tangential 

Velocity 

Fixed 

Free 

Analysis 

X 


X 


Design 

X 

X 


X 


Table I. Differences Between Analysis and Design Regions 


The index PTYPE identifies the type of region. Both 
PTYPE = 0 and 1 designate analysis regions. The difference is 
that a region of PTYPE = 1 is free to translate in the (x,y)- 
coordinate system while satisfying continuity requirements, 
whereas the location of a region of PTYPE = 0 will remain fixed. 
The freedom to translate is sometimes useful if the analysis 
region is nestled between two design regions. Of course, no 
analysis region is allowed to rotate, change shape, stretch, or 
shrink. PTYPE = 2 and 3 designate design regions. A region of 
PTYPE = 3 is free to change geometry while satisfying continuity 
requirements with the adjacent regions. PTYPE = 2 is the same as 
PTYPE = 3, except that for PTYPE = 2, two degrees of freedom are 
eliminated by prescribed fixed values for the (x , y) -coordinates 
of point i A (p). A region of PTYPE = 0 cannot be immediately 
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followed by a region p of PTYPE = 2, because continuity require- 
ments alone would fix the location of i A (p). In that case, 

PTYPE (p) should be selected as 3. Table II defines the allow- 
able ordering of regions on an airfoil element. 


If PTYPE (p) is 

* 

then PTYPE (p+1) can be 

0 (Analysis) 

1 (Analysis) 

2 (Design) 

3 (Design 

0 or 3 only 

1 or 3 only 

0 , 1 , 2 , or 3 
0 , 1 , 2 , or 3 


* If p = p^ , replace (p+1) by 1. 


Table II. Permissible Ordering of Region Types 


There are two additional restrictions. First, at least one 
region per element must be PTYPE = 0 or 2 in order to prevent 
the entire element from unrestricted translation in space. 

Second, the user must divide the element into regions in such 
manner that any sharp corners are located at the intersection 
of two regions. This second restriction simplifies the indexing 
but does not limit the solution generality of the method. That 
is, it is permissible for adjacent regions to be both analysis 
or both design, as identified in Table II. 

It is convenient to designate a set of adjacent regions that 
are free to move as a "group of free regions". On any airfoil 
element q, each group of free regions is assigned an index IGP 
(1 <_ IGP <_ NGP) . The regions included in group IGP are 

P a (!GP), p A +l, ..., P B -1, P B dGP). 


The locations of i^.(PA) and iB^PB) are fixed, but all inter- 
mediate points are free to move. If the trailing edge is in 
the midst of a group of free regions, then the regions included 
in the group are ordered as follows: 

P A' P A +1/ P T -1, p t ' 1 ' 2 ' Pb -1 ' P B- 

Each region of PTYPE =1, 2 , or 3 will be part of a group of 
free regions. The assignment of group indexing IGP, p^(IGP) 
and p B (IGP) is performed automatically by Program MAAD. 
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The transformation of a design region during an iteration 
cycle is illustrated in figure 9. The free stream velocity 
with respect to the (x,y) -system (i.e., and a) will not change. 
Although the normal direction of each panel changes, the value 
remains equal to the prescribed value • The change in 

shape is completely defined by the array 66j_ fi A <_ i <_ (i B ~l) ] . 
The ratio of surface distance of the region after transformation 
to the distance before transformation is (1 + SSF) . It is noted 
that the percentage length change is the same for every panel 
in the region. Therefore, the prescribed tangential velocity 
array V>r ■ describes a fixed distribution of velocity versus 
1 P 

normalized surface distance m the region. 



Figure 9. Transformation of a Design Region Between Successive Iteration Cycles 



The displacement between iterations of point i of design 
region p in figure 9 is 

(i-1) 

6x . = Sx . + E As. { (1 + 6 SF ) cos (0 . + 60.) - cos0 . } (35) 

i i A i=i 1 3D 3 

A 3 i A 

(i-1) 

5y. = 6y. + E As. {(1 + 6SF) sin(9. + 69.) - sin0 . } (36) 

A j=i A J 333 

Equations (35) and (36) can be linearized. 


(i-1) 

6x . = 6x. + E As. {- sin0 . 60 . + cos© . 6SF} 


A j=i, 


3 3 


(37) 


(i-1) 

6 y. = 6y. + E As. {+ cose . 60. + sine. 6SF} 


A 3 = 1 


A 


3 3 


(38) 


By applying the indexing nomenclature for the group of free 
regions associated with design region p, the terms 6xj_ A and 6yj_ A 
in equations (37) and (38) can be expressed as functions of the 
geometry change of preceding regions. Equations (37) and (38) 
become 


P 

6x . = E 
1 

Pi 


E As. {-sin0. 60. + cos0 . 6SF} 


3 = i, 


3 3 


(39) 


6y. = E E AS. {+cos0 . 60. + sin0 . 6SF} 


Pa 3 = 1, 


3 3 


(40) 


where 1q = ( i ,3 - D if the region under summation is not p and 
i^ = (i-1) if the region is p. Equations (39) and (40) apply 
not only for design regions but for analysis regions of PTYPE = 1 
as well. Of course, if PTYPE (p) = 1, the values 60 j and 6SF for 
region p are identically zero and are ignored. 

The complete set of geometric unknowns in a design problem 
is set up as an array 5v j , ordered as follows: 
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6 v 

6v 

6v 


1 

2 

3 


60 . 

69 (i A+ D 


60 


▼ 


(V 1 ) 


6SF 


60 . 

1, 


66 


(i A +1 > 


▼ 


60 


(i B _1) 


6SF 


1 

] 


First region 
of PTYPE 
= 2 or 3 


Second region 
of PTYPE 
= 2 or 3 



Element 1 


Element q 


Element q T 


The method for establishing the first-order relationship 
between velocity and geometry perturbations will now be described. 
The procedure is to repeat every step described above in the 
analysis-only solution formulation. This time, however, each 
quantity is expanded by its differential with respect to the 
variable 6vj c (defined above) . The resulting differential 
analogue of equation (30) is in the form 


Z 

j 


{C. .du . + u . E 
1D 3 3 k 


3 C . . 

^ d V 1 


E 

k 


3 RHS . 



dv k } 


(41) 
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In equation (41), Cij and uj are obtained from the previous 
analysis-only solution. The differentiation of Cij and RHSi is 
performed in three steps. First, derivatives are established with 
respect to £p , np» and s. Second, differentials d£p and drip are 
expressed in terms of d0j_, dx£, dyi, dxi+i, dyi+1, dxp, and dyp 
(where " i" refers to an influencing panel and "P" to a control 
point). Third, equations (39) and (40) are used to express the 
differentials in terms of d0j and dSF, i.e., in terms of dv^. 

Differentials of the induced potential functions (19) - (22) 
are provided in the form 


d ( ) 


d<f> 


d<j> 


d<j> 


51 
VI 

52 


3 ( ) a r . a ( ) ... . a ( ) 

3?p dC P 3n p dn p 3 s 

ds 

(42) 

h lF A d h + F B dri P + F c ds) 


(43) 

b tF B d 5p - F A d "p + [F D - 

(8 + -j) Jds} 

(44) 

b ll ' 1 + & f a + & V d «P 


+ [ «7» f b - i TT ) V 

dn P 


+ [+ < 2? ) ‘ b <F E + 

i ^ p 

4> < 2 > F B ] ds> 

s 

(45) 


d *V2 - b ' 1 <1^ F B - <ir> V d «P 


+ [+ 1 - <r-> f a - ^ f b ) d + 


n P 1 

+ ( 2S ) " ~ (F ir + T) + ( 


SpTlp 


2 } F A ] ds} 


(46) 


From equations (7) and (8) , 

d£ p = tripde^^ + cose^^ [dx p - ^(dx i + dx i+p ) ] 


+ sine i [dy p - j(dy ± + dy i+1 ) ] 


(47) 
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E^de . + cose . 
p i i 


drip = - 


[dy p - jtdy. + dy. +1 )] 


sine . [dx„ 
1 P 


2 (dx i + dx i+ l )] 


( 48 ) 


Equations (39) , (40) , (47) , and (48) are used to express 

dEp and dnp in terms of the independent differentials, d0 . and 
dSF . 3 


Following evaluation of the derivatives, equation (41) is 
expressed in the form 


E C . . du . = z D . . dv, 
j 13 : k Ik k 


(49) 


The system of linear equations (49) is pre-multiplied by 
the matrix inverse of , with the following result: 


dU >! ‘ 1 ([C U ] 


"I 


l <D ik dV k ,) 
k 


(50) 


l (G ik dv k> 

k 


where the array is defined by equation (50) . 

The final desired velocity-geometry perturbation relationships 
are generated by differentiating equations (31) and (32) with 
respect to 0j_, ©ei , Yi, and Y i+i. With the aid of equation (33) 
and the fact that dSF and d6j are equivalent to dvk, plus the fact 
that equation (50) provides the relationship between vortex and 
geometry diffentials, it is possible to express dVTi an< 5 dVTEi i n 
the following form. 

dV m = z A. .dv . ( 51 ) 

T ± j 11 3 ( ’ 


dV T = 1 A E dV i (52) 

j ij J 

The calculated arrays Aj_j and Aei^ are exactly the velocity 
partial derivatives of the panel method analysis solution. For 
example, suppose that the analysis problem has been solved by 
Program MAAD for two geometries, where the difference in geome- 
tries is defined by arbitrary finite values in the array 6vj and 
where the difference in calculated tangential velocity components 
of panel i is 6VTj_* Then the limiting conditions as all values 
6vj approach zero will be 
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<5v .->-0 
3 


E A. . 6v . 

j 1D 3 


In a design problem, values and Vt e . are prescribed 

for each region of PTYPE = 2 or 3. The solution approach is to 
use equations (51) and (52) to find the corresponding geometry- 
perturbation. It may appear that one could simply insert desired 
values 6VTi into the left-hand-side of equation (51) and then 
solve for 6vj. However, Aj_j is usually an ill-conditioned matrix. 
Furthermore, there is no guarantee that geometric continuitv re- 
quirements would be met. For those reasons, it is more appropriate 
to minimize the function E. 


Ss . 

I i { ,1 '5T’ [ Y - Y - l A ij *V 


q p l=l A T 


i 3 


B As . + As . 

+ * (~ - o -5 — ) [V n 


V T ' 2 A E.. 

3 ID J 


2 (i B _1) 

" 2 
+ [ (— ) E As.] 6SF } 

^ r^i 1 » 1 

T l-l A 


Pb ‘V 1 ’ 

+ q IGP YlGP { p !p iS K Sin8 k S6 k + cosS k 5SF)) 


, P B ^B 

+ E E A { E E As, (+ cose, 69, + sine 6SF) } 

q IGP p=p A k=i A k k K k 


Equation (54) requires explanation. The equation expresses 
the mean-square difference between calculated and prescribed 
tangential velocity distributions, with Lagrangian constraints 
for satisfaction of continuity requirements. Only regions of 
PTYPE = 2 or 3 are to be included in the summation. Sip is de- 
fined as the perimeter of element q. In the term 



29 


A s j_ and A si — 1 are to be neglected if i = i B and i A , respectively. 
At sharp corners, VteIp is ignored. A weighting factor W 2 
is selected for each design region by the user to control the 
stretching. A very small value such as 10 -5 has no significant 
effect. However, a large value such as 10 + 5 w in suppress the 
length change. An additional option available to the user is 
the stipulation that two or more successive design regions in a 
group of free regions share a common value for 6 SF . This is 
achieved by consolidating the coefficients of ASF for the appro- 
priate regions and then removing the extra unknown (s) from the 
array Avj . 

The unknowns in equation (5 4) are Avj (which is equivalent 
to <$e k and ASF in accordance with the definition of 6 v j ) and the 
Lagrangian multipliers A X/ jqp and ^y,lGP* The unknowns are 
determined by setting equal to zero the derivative of E with 
respect to each unknown and then solving the resulting system of 
linear equations by a standard matrix method. 

After the array Avj is calculated, the geometry is corrected 
one group of free regions at a time. Starting with panel i A of 
region p A , the panel is stretched and rotated the required 
amount. This provides new coordinates for point i A +l. The pro- 
cedure is repeated for panel i A +l, etc., until all panels in the 
group of free regions are included. Then the new coordinates 
of point ip of region pp are examined to verify that continuity 
requirements are satisfied. Even though the Lagrangian con- 
straints of equation (54) provide satisfaction of continuity 
requirements, the constraints are only first-order accurate with 
respect to Avj . Define the discontinuity vector (DXB , DYB) as 


: . , . ] - [x . . . ] 

1 B P B NEW 1 B iP B J REQUIRED 


'B ' NEW 


_ [Y i (P ) ] 

B vp B ; REQUIRED 


If the magnitude of ^DXB^ + DYB 2 is greater than some very small 
tolerance, it is necessary to make a closure correction to the 
geometry. What is desired is a new set of values A 9j_ and ASF 
for each design region that will bring ^DXB 2 + DYB 2 to within the 
prescribed tolerance. It is important for the closure correction 
to change the geometry as little as possible. This is accomplished 
by minimizing a new error function E. 
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P R (i R ~D 

B B 2 2 2 

E 5 Z {[ Z AS. (60.) ] + (W +1) (DST ) ( 6SF) } 


p=p a 1=1, 


1 X 


B 


(V 1 ) 


+ X { DXB + Z [ (DX ) 6SF - Z As. sin9 . 60.]} 


P=P, 


1=1 


A 


1 11 


B 




+ A {DYB + Z [ (DY ) 6 SF + Z As i COS0 ^ 66^} 


P=P 7 


1 = 1 , 


(57) 


where 


DST 


(i B (p)-D 


i = i A (p) 


A s . 
1 


(58) 


DX 


(i B (p) -1) 

z 

i=i A (p) 


a s . cose . 

i i 


(59) 


DY 


(i B (p)-l) 


i=i A (p> 


A s . sine . 
l l 


(60) 


In equation (57) , only regions of PTYPE = 2 or 3 are to be in- 
cluded in the summation over p. The unknowns are 6 0-^, 6SF, A x , 
and Ay/ where A x and Ay are Lagrangian multipliers. The first- 
order Lagrangian constraints correspond to the elimination of the 
discontinuity vector (DXB, DYB). The weighting factor W 2 and the 
option in which successive design regions share a common length 
variation 6SF apply to equation (57) in the same manner described 
for equation (54). 

E is to be minimized with respect to the unknowns. 

0 = e ^ q- \ = (2As . ) 60 . - (As. sine.) A + (As. cose.) A (61) 
3 ( o 0 j ; 1 j 3 jx i jy 

0 = TtUfT = [ 2 -DSTp. (W 2 + 1)] 6SF + DX p -A x + DY p -X y (62) 
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A mathematical artifice which will result in reduced computing 
effort is to multiply equations (61) and (62) by appropriate fac- 
tors and then sum the products. 


P B • ‘V 1 ’ 
0 = Z { [ E 

P=P A j=i A 


sine . 
3 


E 

(60 j) 


+ 


DX 

E 

DST (W 2 + 1) 
P 


9 E 


9 (6SF) 


(63) 


0 = 


B 

Z 

P=P 7 


( [ 


Z 

j = ± A 


+ cose 


9 E 


9(66.) 

D 


DY 


] + 


DST (W 
P 


+ 1) 


3 E 

3 ( 6 SF ) 


} (64) 


By substituting the right-hand-sides of equations (61) and (62) 
into equations (63) and (64), and then by applying the fact that 
the Lagrangian constraint functions of equation (57) are equal 
to zero, equations (63) and (64) become 


2 • DXB = 


B 

+ X Z 
x 

P=P Z 


(D V 


DST (W + 1) 
P 


(i B _1) 
+ Z 


2 

As . sin 0 . ] 
3 3 


+ X Z 

Y 


p=p 


B (DX ) (DY ) 

- E E 

A ‘DST p (W 2 + 1) 


(i B _1) 


3 = 1 7 


As. sine, cose.] (65) 
3 3 3 


2 • DYB = 


r B 

+ X Z 

x 

P=P 7 


(DX ) (DY ) 

- E E 

2 

DST (W + 1) 

P 


(i B -1) 

Z 

3 =i A 


A s . sine . cose . J 
3 3 3 


B 

+ X Z 

Y P^P, 


(DY p } 


DST (W + 1) 
P 


(i B- 1) 

Z 

3=i^ 


As . COS 0 . ] 
3 3 


( 66 ) 


Equations (65) and (66) provide a 2 x 2 system of linear equa- 
tions which are solved for X x and X y . Then equations (61) and 
(62) provide values immediately for each of the other unknowns, 
66 j and 6SF. 
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Using the set of values for 69 j and 6SF, the geometry is 
corrected and a new discontinuity vector (DXB,DYB) is generated. 
The closure correcti on procedure is repeated as many times as 
necessary to bring ^DXBZ + DYB2 to within the prescribed toler- 
ance. Usually, one or two times are sufficient. 

It is noteworthy that the closure correction procedure is 
also applied to the starting geometry in design problems before 
any aerodynamic calculations are performed. This significantly 
simplifies user input requirements. As illustrated in figure 1, 
large discontinuities are permitted. Program MAAD will auto- 
matically generate the desired closure without changing the 
relative panel length distribution in each design region. 

Following the closure correction in a design problem iter- 
ation cycle, a test is performed to determine whether any two 
panels intersect improperly. The only acceptable intersections 
are, of course, between panels i and i-1 at point i. Improper 
intersections destroy the separation between the real and imagi- 
nary flow fields, thereby violating the theoretical requirements. 
There can be situations in which the user will have to redefine 
either the starting geometry, prescribed tangential velocity 
distribution, or gap between airfoil elements. However, Program 
MAAD has an option for increasing the thickness of an airfoil 
element if the upper and lower surfaces overlap. This option is 
usually sufficient to eliminate improper panel intersections. 

At this stage, the inverse step of a design problem itera- 
tion cycle is complete. All data necessary for the analysis step 
of the next iteration cycle has been calculated. 

Prediction Accuracy of Analysis Solutions 

In references 1 and 2, predictions by the analysis method of 
Program MAAD (using constant source panels) were compared with 
predictions by both a low order and a higher order method that 
impose flow tangency boundary conditions directly. It was con- 
cluded that the MAAD method is competitive with the higher order 
formulation for a wide variety of geometries. Rather than repeat 
the detailed analysis examples, a brief discussion on the rela- 
tive prediction accuracies of the three methods is presented 
here . 


Figures of merit for prediction accuracy are the calculated 
drag and the difference between lift calculated by pressure 
integration and circulation for solid, closed airfoils. Analy- 
tically, both the drag and lift difference are zero. For single 
element airfoils modelled by from thirty to forty panels, the 
drag coefficient is typically on the order of 0.0005 for both 
the MAAD method and higher order method, versus 0.0025 for the 
low order flow tangency method. Similarly, the lift difference 
is typically 0.001 for the MAAD and higher order methods, versus 
0.005 for the low order flow tangency method. 
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For a circular cylinder represented by twenty equal length 
panels, the surface distribution of relative error in calculated 
midpoint velocities is essentially uniform for each of the three 
methods. Specifically, the error is 2.06% for the low order flow 
tangency method, 0.22% for the MAAD method, and 0.05% for the 
higher order method. 

It is noteworthy that there is seldom any appreciable 
difference in prediction accuracy between the use of constant 
and linear source panels in Program MAAD . Furthermore, the 
extrapolated bisector Kutta condition (using the recommended 
value k 2 = 0-02) will generate a lift coefficient that is typi- 
cally within a tolerance AC^ < 0.01 of the mathematically exact 
value, whether the panels in the MAAD method have constant or 
linear source distributions. This 0.01 tolerance is based on the 
use of from 30 to 40 panels to model a single element airfoil. 

The reason behind the good prediction accuracy of the MAAD 
method centers on the use of internal potential boundary condi- 
tions. A flat panel source-vortex model characteristically under- 
predicts tangential velocity at panel centers, but overpredicts 
near endpoints. On the other hand, the average velocity between 
adjacent panel centers, which is the difference in potential, is 
significantly more accurate. 

Although it is true that the MAAD method exhibits higher 
order prediction accuracy in most practical problems, there are 
exceptions. For example, the predicted tangential velocities 
induced by positive blowing on both the upper and lower surfaces 
of a thin, cambered airfoil will exhibit low order accuracy . 
However, it is believed that such limitations are outweighed 
by the simplicity of the basic formulation, especially when 
applied to the more complex problem of design. 

Example Design Solutions 

Five design solutions are presented to demonstrate the 
flexibility of Program MAAD. It is noted that in the example 
solutions, constant source panels were used in the program. 

The linear source feature described in the preceding formulation 
sections was added after the figures were prepared. However, 
some of the examples have since been repeated using linear sources, 
and there are no significant differences in the results. 

A typical design solution requires five iteration cycles, 
where each cycle requires 100%-300% more computing time than an 
analysis-only solution. On the CDC CYBER 175, a two-element 
airfoil represented by a total of 70 panels requires 20 seconds 
computing time per iteration cycle for a design-only solution. 

The program can handle as many as five airfoil elements simul- 
taneously . 
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The selection of examples is intended to demonstrate design 
solution characteristics for a variety of applications . Included 
are two design-only solutions for two-element airfoils. The 
design of a cavitation body illustrates the value of the variable 
panel length feature. The ability to modify airfoils to offset 
the effects of wind tunnel interference is demonstrated for a 
Liebeck high-lift geometry. A simple separated flow model is 
presented in which the pressure gradient in the recompression 
region is calculated by the program. 

The first example was solved previously by the analysis-only/ 
design-only method of references 1 and 2. The identical solution 
by Program MAAD is presented here to demonstrate prediction 
accuracy and convergence properties. In figure 10, the analysis- 
only solution for a two-element airfoil is compared to an exact 
analytical conformal mapping solution by Williams (reference 8) . 
Notice that the lift predicted using the extrapolated bisector 
Kutta condition ( K 2 = 0.02) is quite accurate. Not shown in 
figure 10 is a comparison made at the exact Ci = 2.03 between 
calculations by Program MAAD and the higher order analysis method 
of reference 5. Over the complete two-element airfoil boundary 
modelled by 66 panels, the average magnitude of the difference 
in calculated pressure coefficients is less than 0.005, demon- 
strating the good agreement between Program MAAD and the higher 
order formulation. A design-only problem was established by 
unwrapping both airfoil elements to provide the open rectangle 
starting geometry of figure 11. The velocity distribution 
calculated by Program MAAD for the analysis solution of figure 10 
was prescribed on the starting geometry. A very large length 
weighting factor was selected for each region to prevent any 
change in panel lengths. The converged geometry generated by 
Program MAAD in five iteration cycles is virtually identical to 
the target geometry of figure 10. Specifically, the coordinates 
at all 66 panel endpoints differ by less than one thousandth of 
one percent of the main element chord. 

In the second example, the problem is to design a symmetric 
solid airfoil with a flat plate nose and tail such that at zero 
incidence the upper and lower surfaces have a constant velocity 
V = 1.50 V rj (Cp = -1.25) . A square was selected as the starting 
geometry (figure 12). The flat plate nose and tail are analysis 
regions, selected as PTYPE = 0 and 1, respectively. The upper 
and lower surface design regions are both PTYPE =3. A very 
small weighting factor (W = 10"^) was selected for both of the 
design regions in order to allow unrestricted stretching. The 
design solution converged to within plotting tolerances in five 
iteration cycles, and the calculated velocities agree well with 
the prescribed distribution. The ability of Program MAAD to 
calculate both shape and surface distance was essential to the 
design of this cavitation geometry, for which the chord was unknown. 

In the third example, the geometry of a Liebeck airfoil was 
modified to offset non-uniform flow angularity induced by simu- 
lated wind tunnel walls. The modification was conducted for a 
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Figure 10. Two-Element Airfoil Analysis Solution 





(Indistinguishable from Target Geometry) 


Figure 1 1 . Two-Element Airfoil Design Solution 
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Converged Solution 



Figure 12. Design of a Constant Pressure Geometry 


test section height- to-airf oil chord ratio of two. Although much 
larger ratios are standard in practice, here the objective is to 
evaluate the design method for very strong interference. The 
airfoil was positioned approximately midway between the walls and 
pivoted at the upper surface quarter chord point. Blockage ef- 
fects were accounted for by applying an average channel velocity 
equal to 0.9822 V re f, where V re f is the reference velocity for 
Cg and Cj,. The three pressure distributions presented in figure 
13 were calculated by Program MAAD at the common lift coefficient 
1.80 , corresponding to the Liebeck airfoil in unbounded flow 

at 11.2° incidence (reference 9). In spite of correcting the 
channel velocity for blockage and correcting the chord incidence 
for the average induced flow angularity, the pressure distribution 
of the unmodified Liebeck airfoil in the wind tunnel reveals signi- 
ficant wall effects. The design solution was generated by pre- 
scribing the unbounded flow pressure distribution and then solving 
for the geometry (including chord incidence) in the presence of 
tunnel walls with channel velocity 0.9822 V re f . To within plotting 
tolerances, the design solution completely offsets the effect of 
tunnel walls on the calculated pressure distribution. 

In obtaining the design solution of figure 13, the airfoil 
element was divided into four design regions. Region p=l spanned 
most of the lower surface, region p=2 encompassed seven panels 
surrounding the stagnation point, region p=3 terminated at the upper 
surface quarter chord pivot point, and region p=4 covered the 
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remainder of the upper surface. The length weighting factors 
selected for the four regions were, respectively, W = 10 “^, io - 5, 
10 +5 , 10 + 5. The low weighting on the first two regions allowed 
the arc length location of the stagnation point and the velocity 
gradient near the stagnation point freedom to adjust. This 
approach can be extremely effective in converting what would 
otherwise be a nonexisting prescribed velocity distribution into 
approximately the closest possible existing distribution, without 
disturbing the prescribed velocities. 
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Figure 1 3. Design to Eliminate Wind Tunnel Wall Effects 



The two element airfoil of figure 14 was designed by Program 
MAAD on an interactive graphics system by applying the high 
lift rules of Liebeck and Smith (references 9-11) . Both the 
adaptation of Program MAAD for interactive graphics and the two- 
element design were conducted as part of the 1978 MCAIR Indepen- 
dent Research and Development Program. The design pressure 
distribution was established for a Reynolds number of approxi- 
mately 5 million, C& = 4.7, and a = 25°. The reference length 
is the extended chord C of figure 14. Both elements are designed 
for the steepest possible upper surface pressure gradient without 
inducing turbulent separation (reference 12). This allows for a 
very low pressure upper surface roof-top acceleration region. 

The presence of the short chord rear element raises the trailing 
edge (dumping) velocity of the main element. This allows for 
the magnification of the upper surface pressure distribution 
of the main element beyond that of the rear element without in- 
creasing the likelihood of separation. The lower surface of 
each element was designed to be as close to the maximum stagna- 
tion pressure as possible while maintaining reasonable thickness 
for structural purposes. 
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Figure 14. Automated High Lift Design Solution 
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A simple model for predicting the pressure distribution on 
an airfoil with trailing edge separation is presented in figure 
15. The model is based on semi-empirical rules for the behavior 
of pressure distributions in separated flow. By using the mixed 
analysis-design mode of Program MAAD, the corresponding geometry 
of the wake can be generated. The combined airfoil-wake geometry 
is represented as a closed streamline. The wake consists of a 
constant pressure upper surface region (p=3) followed by upper 
and lower recovery regions (p=4 and 1) . All the regions are 
PTYPE = 3, except for the solid airfoil region p = 2, which is 
PTYPE =0. A very large length weighting factor (W = 10 + ^) is 
selected for region p = 3 . Therefore, the upper recovery region 
p = 4 will begin approximately directly above the trailing edge 
of the actual solid airfoil. The prescribed velocity recovery in 
both regions p = 1 and 4 varies linearly from the separation 
velocity level to a nominal value slightly less than free 
stream. This will generate zero wake loading. A very small 
length weighting factor (W = 10“^) is selected for regions p = 1 
and 4. Consequently, the wake length is free to change. The 
program will indirectly generate the recovery pressure gradient 
by determining the coordinates of the wake trailing edge that are 
consistent with the entire model. This is in contrast to the wake 
model established by Henderson (reference 13) , in which the 
recovery pressure gradient is essentially fixed a priori. 
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Figure 15. Wake Model for External Inviscid Flow Field 
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The application of the wake model of figure 15 to a NACA 
633-018 airfoil at 15.6° incidence is illustrated in figure 16. 

The wind tunnel data of reference 14 was used to establish both 
the separation point location and wake velocity V w for the wake 
model. The agreement between calculated and experimental 
pressures is good, indicating that the present model is reason- 
able. It is noted that the effects of the boundary layer dis- 
placement and the wake mass flux deficit were neglected. It is 
anticipated that by coupling the equations of integral viscous 
flow theory with the existing inviscid velocity-geometry pertur- 
bation functions of Program MAAD, it will be possible to accurately 
solve for strong viscous-inviscid interactions without resorting 
to empirical rules. 
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Recommendations and User Guidelines 

Program MAAD is an accurate and reliable method for handling 
a wide variety of two-dimensional, incompressible, potential flow 
problems. Basic user guidelines and recommendations for further 
program development are briefly described here. 

Design solution convergence is usually very good, provided 
that the prescribed velocity distribution is reasonably close to 
an existing distribution. For every geometry there is an exist- 
ing pressure distribution; however, the converse is not true. 
Fortunately, in the iteration process, the program will generate 
a velocity distribution in the analysis step that will guide the 
user toward necessary modifications in the prescribed distri- 
bution. This guidance is especially helpful if interactive 
graphics are used. 

In design-only problems, the user is often faced with the 
dilemma of where to prescribe the stagnation point and how steep 
to make the local velocity gradient. There are two approaches. 
First, let the program run through two or three iteration cycles. 

If the stagnation location and local gradient are inconsistent with 
the remainder of the prescribed distribution, the designed leading 
edge region will appear saw-toothed. The user can then smooth 
the leading edge and run an analysis solution to determine an 
existing distribution to prescribe in the vicinity of the stag- 
nation point. The second and simpler approach is to divide the 
element into three or more regions, one spanning the lower sur- 
face with low length weighting, one in the vicinity of the stag- 
nation point with low weighting, and one covering the upper 
surface with high weighting. As described in the third solution 
example, this second approach will automatically generate the 
proper adjustment to the stagnation region. 

In the modelling of bounded flow problems such as wind 
tunnel simulation, there are two rules. First, the net velocity 
flux through all boundaries must be zero, which is a direct 
consequence of the continuity equation. Second, it is desirable 
to define values for V m and a that approximately correspond to the 
average anticipated fluid velocity. Theoretically, any values are 
acceptable because there is no physical free stream in a completely 
bounded flow.. However, prediction accuracy and design solution con- 
vergence are usually enhanced if reasonable values are selected. 

It is recommended that Program MAAD be operated on a computer 
system providing twelve or more significant figures. This is 
important only for exceptionally thin airfoils. It can be shown 
mathematically that the calculated loading is a function of the 
difference in potential between adjacent upper and lower control 
points. If the control points are sufficiently close, then the 
difference in potential can be several significant figures less 
accurate than the potential. For this reason, the possible 
future introduction to the program of intermediate and far field 
influence coefficient expansions should be approached with 
caution. 
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Although there are many conceivable improvements that could 
be made to Program MAAD with respect to both extended capability 
and computing efficiency, two are most worthy of mention. The 
first is to extend the method to handle multi-energy flows such 
as jets. An obvious approach would be to solve for the combined 
source-vortex distribution and geometry that simultaneously 
satisfy internal and external boundary conditions for the two 
(or more) flows of different potential. Another approach that 
would not increase computing effort significantly but would re- 
tain the benefits of Green's third identity would be to handle 
a real and an imaginary internal flow in concert with a real and 

an imaginary external flow. The second conceivable improvement to 
Program MAAD is to implement an accurate method for handling 
airfoils with finite thickness trailing edges, which is the sub- 
ject of the next section. 

KUTTA CONDITION FOR FINITE TRAILING EDGE THICKNESS 

For airfoils with sharp trailing edges, the previously des- 
cribed extrapolated bisector Kutta condition of the Multi-Element 
Airfoil Inviscid Analysis and Design Program (MAAD) typically 
provides calculated lift coefficients that differ by less than 
0.01 from the exact analytic value. Even this small error will 
vanish as the number of panels is increased indefinitely. How- 
ever, many advanced airfoils have trailing edge thicknesses of 
1/4 to 3/4 percent chord, and for these airfoils Program MAAD 
characteristically overpredicts the lift by ten to twenty per- 
cent. Furthermore, the error does not diminish with increased 
panel density. The overprediction results from the common prac- 
tice of not panelling the base. Consequently, the local flow is 
not properly modelled and the correct circulation is not gen- 
erated. A method for modelling the base has been established 
which typically improves the lift predictions by a factor of ten. 

The objective of the method is to use base panelling to 
generate smooth flow off both the upper and lower surface trail- 
ing edges. The extrapolated bisector Kutta condition is replaced 
by the base panel control equations defined in figure 17. The 
upper and lower surface velocities immediately upstream of the 
trailing edge are Vg and Vg, respectively. The velocities will 
be tangent to the surface unless nonzero normal velocity boun- 
dary conditions are prescribed. The velocity component in the 
direction of the trailing edge bisector (Vg) is to vary linearly 
from Vbu at the upper surface to Vbl at the lower surface. Here 
the trailing edge bisector is defined as the line intersecting 
the mid-point of the base at an angle equal to the average of the 
upper and lower surface trailing edge angles. 

The number of panels used to model the base is not especi- 
ally significant. One or two is usually adequate. It is 
recommended that the values for Vg , Vg, Vgg, and V B g be calcu- 
lated by linearly extrapolating the mid-point velocities of 
neighboring panels to the upper and lower surface trailing edges. 
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Figure 17. Base Panel Control Equations 

The base panel control equations have not been coded expli- 
citly for Program MAAD. Therefore, it is presently necessary to 
satisfy them manually by the proper combination of linearly 
independent velocity distributions for a unit free stream, unit 
circulation, unit constant base blowing, and unit linear base 
blowing. The proper combination is uniquely determined by the 
base panel equations of figure 17. This approach was used in the 
examples presented below. 

The examples compare the inviscid solution predictions of 
four methods applied to airfoils with finite thickness trailing 
edges. The methods are (1) Catherall-Sells (reference 15), (2) 

Oellers (reference 16), (3) Program MAAD with the extrapolated 

bisector Kutta Condition and no base panels, and (4) Program MAAD 
with the present base panel modelling. The Catherall-Sells 
conformal mapping method is virtually exact for sharp trailing 
edge airfoils and provides reasonably smooth flow at the corners 
of a blunt base. Therefore, Catherall-Sells is considered to 
be a standard of comparison. The Oellers vortex-only panel 
method does not provide smooth flow at a blunt base, but the 
Oeller equal pressure upper surface-lower surface Kutta condi- 
tion generates reasonably accurate lift predictions. The 
geometries examined were the three simple symmetric airfoils of 
figure 18 and the Whitcomb advanced airfoil of figure 19. 
Identical input coordinate definition was used for each of the 
four prediction methods. The base was left unpanelled except 
for the Program MAAD calculations using the present base model- 
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ling. In every example, piecewise constant source distributions 
were used in Program MAAD. In each of the four methods, the 
lift coefficient is referenced by the (maximum) airfoil chord 
and was calculated by integrating the pressure force acting 
normal to V^. 
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Figure 18. Symmetric Airfoils Selected for Blunt Base Analysis 





t XE = 0.0071 C 


Figure 19. GAW-1 Airfoil 
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For each of the three symmetric airfoils of figure 18, 120 
panels were used to represent the geometry forward of the 
trailing edge. The panelling is sufficiently dense that the 
difference in the predictions of the four methods should be 
dominated by the base modelling and not by numerical discretiza- 
tion. Depending upon the trailing edge thickness, the base was 
represented by either 10, 5, or 0 panels. Solid body boun- 
dary conditions were prescribed ahead of the base. 

Plotted as a function of base (trailing edge) thickness in 
figure 20, the predicted lifts for all four methods attain a 
common value as trailing edge thickness vanishes. However, the 
lift is significantly over-predicted by Program MAAD in the 
cases for which a finite base was left unpanelled. Clearly, the 
present base modelling alleviates the problem. The accuracy 
is also evident in the pressure and velocity distributions of 
figures 21 and 22. 



Base Thickness - tyg/C 


Figure 20. Lift Prediction for Symmetric Blunt Base Airfoils 
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Figure 22. Blunt Base Airfoil Analysis (90° Incidence) 

For the GAW-1 Airfoil of figure 19, lift prediction error is 
plotted against angle of attack in figure 23. Here the error is 
defined as the difference between the calculated lift and the 
lift predicted by Catherall-Sells . With only two panels, the 
base modelling reduces the Program MAAD error by an approximate 
factor of ten. 

Although the present base modelling method was developed for 
Program MAAD, the method can be implemented directly in any 
surface panelling program. The modelling is a simple but effec- 
tive approximation to the fact that no real fluid can turn a 
sharp corner on an impermeable surface. 

THREE-DIMENSIONAL POTENTIAL FLOW METHOD 

In 1978 a surface panel-method computer program was devel- 
oped for solving incompressible potential flow analysis problems 
for arbitrary lifting or nonlifting three-dimensional geometries 
(reference 1). When applied to thin wings, the program demon- 
strated significantly better prediction accuracy than the Hess 
(Douglas Neumann) source method (reference 17) . This improve- 
ment was the anticipated result of using the mild source-doublet 
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Figure 23. Lift Prediction Error for GAW-1 Airfoil 
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distribution of Green's third identity in place of a distribution 
having strong source density. However, the prediction accuracy of 
the program was dependent upon the degree to which neighboring can- 
els were uniform in size and shape. Furthermore, the modeling was 
limited to trapezoidal panels. These restrictions were elimin- 
ated as part of this contract by modifying and extending the 
earlier program, thereby creating the MCA1R 3-D Subsonic Poten- 
tial Flow Program (Version 1) , hereinafter designated the 
"MCAIR 3-D Program" . 

As in the earlier program, the geometry in the MCAIR 3-D 
Program is modeled by flat, quadrilateral panels with distri- 
butions of constant source and quadratic doublet density. 

Potential influence functions from reference 18 have been 
incorporated to eliminate the trapezoidal recuirement 
of the earlier functions. The quadratic doublet distribu- 
tion on each panel is established by a least squares spline 
fit through neighboring panel control points. In the earlier 
program, only nine control points were used to establish the 
spline fit, resulting in doublet discontinuities at panel edges 
and subsequent solution sensitivity to nonuniform paneling. In 
the MCAIR 3-D Program, twenty-one control points are used in 
such a manner that the spline is nearly continuous at panel 
edges. This spline is similar to the one developed for the 
Boeing PANAIR Program (unpublished); however, in the present 
case, the spline can be carried across the intersection of net- 
work edges, even if the intersection includes a wake. This 
eliminates the requirement for edge boundary condition control 
points on smooth geometries such as fuselages and guarantees 
that axisymmetric flow problems will have axisymmetric solutions. 
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Presented in this section is a description of the solution 
method. Experience with the analysis mode of the Multi-Element 
Airfoil Inviscid Analysis and Design Program guided the develop- 
ment of the MCAIR 3-D Program. In that regard, the present 
three-dimensional analysis method can be considered to be a 
nearly direct extension of the two-dimensional formulation 
described in a preceding section. In accordance with Green's 
third identity, the source density on each panel is determined 
by the orientation of the panel with respect to the free stream 
direction. The doublet distribution is calculated by requiring 
that the internal perturbation potential be zero at selected 
control points, corresponding to a uniform imaginary internal 
flow. It is recommended that the reader be acquainted with 
the theoretical and numerical background presented in refer- 
ence 1. The formulation for the MCAIR 3-D Program is described 
below, to be followed by a simple wing example solution. 

Geometry Panel Modeling 

The surface of each body in the flowfield is divided by 
the user into one or more sections (figure 24) , where each sec- 
tion will be further subdivided into panels. A typical section 
is a wing, tail, or part of a fuselage. A body is defined as 
any closed surface with an interior, such as an airplane, 
sphere, or toroid. Canard-wing-tail geometries and slotted 
flap wings are permitted. However, the user must define the 
wake trajectory a priori. The division of a body into sections 
must be performed such that one side of each section is wetted 
by the real external flow and the other by the imaginary internal 
flow. If there are one or two planes of symmetry, only the 
geometry on one side of the plane (s) is defined, and the image 
geometry will be generated automatically. 

With the exception of the following two restrictions, the 
division of a body into sections is arbitrary. First, a section 
cannot be multiply connected in the sense that there are holes in 
the surface. Second, every sharp edge on a body must be aligned 
with section edges, where a sharp edge is defined as the locus 
of points at which there is a discontinuity in the surface 
normal direction. 

Although open wing tips and open fuselage bases are per- 
mitted in the MCAIR 3-D Program, it is recommended that they be 
closed in accordance with the definition of a body. The closing 
seals the boundary between the real external and the imaginary 
internal flows. Of course, physical openings such as flow- 
through nacelles should not be closed. 

The difference between a lifting and nonlifting body is 
dictated by the presence or absence of an intersecting wake. 

A wake is divided into sections in the same sense as a body, 
except that a wake is a sheet and not a closed geometry with an 
interior. No part of a wake can lie in the interior of a body. 


50 




Figure 24. A Section of a Body Surface 

One edge of a wake section will be assigned control points to 
determine the lift. That edge must coincide with a sharp edge 
of a body. Typically, this will correspond to the intersection 
of a wake with a wing trailing edge. The opposite wake edge, 
usually designated the starting vortex, can also intersect a 
body section at a sharp edge. In most cases, however, the start- 
ing vortex will be positioned several body lengths downstream 
to simulate a semi-infinite trailing vortex sheet. Either of 
the two remaining wake section edges can intersect an edge of 
another wake section, intersect a body section, or be free of 
any intersections. The intersection with a body section can 
occur along any path of panel edges selected and identified by 
the user. On an airplane, the inboard wake edge should follow 
the contour of the fuselage, even if the fuselage is tapered. 

Each section is assigned an index i (1 < i < n ) . The 

O o O 

type of section is identified by STYPE(ig) , equal to 1 for a 
body section and 2 for a wake. A third type, designated a 
vortex loop, is also available. The difference between a vortex 
loop and a wake section is that for the former, a single edge 
control point is used to determine a uniform doublet distribu- 
tion. 


The subdivision of a section into panels is performed using 
the rectangular indexing nomenclature of figure 25. It is noted 
that only the indexing is rectangular, not the section or panel 
geometries. If the section in figure 25 is on a body, the 
exterior is assumed to be on the viewer's side of the figure. 
Cumulative point and panel indices are defined as: 
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Figure 25. Indexing Nomenclature for a Section of Panels 


k = k o (i g ) + ( i — 1 ) *n + j (67) 

k p = kp Q (i s ) + (ip-1) • (n-1) + j p (68) 

where k (i c ) and kp (i ) are respectively equal to the total 

number of corner points and panels on preceding sections 1 
through (is-1) . Coordinates (x,y,z) are input for each corner 
point k. With the exception of points on a section edge, each 
point k is a corner of four surrounding panels. At the inter- 
section of two sections, there is no requirement that the sections 
share the same point distribution along the common edge. There- 
fore, it is permissable for a section with dense paneling to 
intersect a section with sparse paneling. Except for sharp 
edges on a body, it is not necessary for the input point distri- 
butions of intersecting edges to lie on a common path of connect- 
ed straight line segments. At sharp edges, boundary condition 
control points will be inserted near the intersection, and there- 
fore gaps or overlaps perpendicular to the intersection could 
adversely affect the calculations. 

The four edges of a section are ordered counterclockwise by 
index i gE , where i gE = {1, 2, 3, 4} respectively corresponds to 
{ j = 1 , i=m, j = n, and i = 1}. Each edge of each section is 
categorized by three values selected by the user, defined below. 
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ETYPA(i gE ) 


0 This edge is either (1) a smooth edge on a body (no surface 
slope discontinuity) or (2) edge number ig E =2 , 3, or 4 

on a wake. In either case, no edge control points will be 
used . 

1 Edge control points will be used. This can be a sharp edge 
on a body section or edge number i SE = 1 on a wake section. 

In the latter case, this edge must intersect two body section 
edges of ETYPA = 1 . 

ETYPB ( i gE > 


0 This edge has finite length. 

1 This edge is a point. This is the only place where triang- 
ular panels are allowed. 

ETYPC (i gE ) 


0 This edge does not lie on a symmetry plane. 

1 The x- z plane is a symmetry plane, with which this edge is 
coincident . 

2 The x-y plane is a symmetry plane, with which this edge is 
coincident . 

A schematic of boundary condition control points on a body 
section is presented in figure 26. There is one control point 
at each panel center. The edge control points are slightly 
retracted toward the center of the section to prevent contact 
with singular edge properties. All body section control points 
are considered to lie on the internal side of a panel. 

For a wake section, all control points are omitted except 
for those on edge igg = 1. Each wake control point is double- 
valued, i.e., there is one point slightly above and one slightly 
below the panel plane. As a Kutta condition, the average velo- 
city between the two points will be set to zero. The direction 
of the vector connecting the two points is input by the user. 

Each panel (i p , j p ) of section i g is described by the four 
corner points (i,j) = (ip, jp) , (i p +l, j p ) , (i p +l, j p +D » 

(i p , j p +l) . In figure 27, these four points have been locally 

reindexed as i = 1, 2, 3, 4, respectively. The edge midpoints 

are i = 5, 6, 7, 8 and the panel center is i = 9. Here the four 
corner points will be adjusted to form a planar quadrilateral 
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ETYPA (3) =£ 1 



Edge Control Points are Used iff ETYPA = 1 


Figure 26. Control Point Locations on a Section 


panel, and a local (£,n,C) - Cartesian coordinate system will 
be defined in relation to the general system (x, y, z) of 
figure 24. 

Using the local indexing nomenclature of figure 27, the 
panel center is determined. 



► 

i - Direction 


Figure 27. Local Indexing Nomenclature and Coordinate System for a Panel 
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ETYPA (2) A 1 




X 9 j 

X 1 + X 2 + X 3 + X 4 1 

[ 1 

y 9 = 4 

Yl + y 2 + y 3 + y 4 

I 

z, + z„ + z_. + z . J 

9 ' 

1 2 3 4 ' 


Define 



as follows: 



p ii 


X i " X 9 J 

-y 



| I 

Hi 

•H 

Oh 

P i2 

III 

| yp " y 9| 

i 

P i3 

1 

Z i - Z 9 * 


(1 < i < 9) 


( 69 ) 


(70) 


Vector is the displacement of point i from the panel center. 


Unit vectors for the (£, n ,?) -system are defined, in conjunc- 


tion with a rotation matrix [a^]. 


e 


(a 21' a 22 ' a 23* 


1P 7 - P 5> 
I P 7 - P 5 I 


(71) 


( a 31 ' a 3 2 ' a 33 } 


(P 6 ~ V * % 

(5 6 - P 8 } X % 


(72) 





(a ll' a i 2 ' a l3 ) 


(e n x a t ) 


(73) 


Unit vector e ? is normal to the panel, directed toward the real 
flow if the panel is part of a body section. The rotation matrix 
[aij] provides the conversion between the local and general coord- 
inate systems. 




] 



(74) 


The coordinates of the nine points in figure 27 are established 
in the (E , , n , C) -system. 
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a lk 




/ 3 





n i H s 

1 k=l 

a 2k 

P ik 

(1 < i < 9) 

(75) 

c i 

' 0 





It is noted that = 0 for all nine points, thereby providing 
the desired planar adjustment. 

The planar adjustment affects the four corner points, 
i = {1, 2, 3, 4}. Consequently, small gaps or overlaps can 
occur between panels. These geometric discontinuities are of 
little significance unless control points are in the immediate 
vicinity, which is the case with edge control points. There- 
fore, the planar adjustment is modified such that a panel edge 
will not be disturbed if that edge is either a sharp edge of a 
body section or an edge of a wake section. The preceding state- 
ment can apply to no more than one edge or two intersecting edges 
per panel . 

Singularity Distributions on a Panel 

The effect of each panel on the flow is represented by a 
constant source distribution and a quadratic doublet distribu- 
tion. The source strength is identically zero on all wake panels. 
For body panels, the strength is determined immediately by the 
application of Green's third identity, 

o = V Np - (76) 

where Vnp is the prescribed normal velocity boundary condition 
at the panel center. The doublet distribution is established by 
fitting a least-square-error quadratic to the doublet strengths 
at the nine points of figure 27. The least squares fit is estab- 
lished before boundary conditions are satisfied, which means 
that values for the doublet strengths must be treated as un- 
knowns . 

The quadratic doublet distribution on a panel is 

u(5,n) = 3-^ + + S 3 n + $ 4 £ 2 + 3 5 £n + (77) 

where the coefficients g-j_ (1 <_ i <_ 6) are constants. If yj 
designates the doublet strength at point j of figure 27 
(1 < j < 9) , it is desired to calculate the matrix b. . such that 

— J — i] 

9 

3, = I b. . y. (1 1 i 1 6) (78) 

j=l 1J J 
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Matrix b^j dictates the influence of each doublet strength y_j 
on each quadratic coefficient 8 . By using the method of least 
squares, bpj will be calculated as a function of the coordinates 
of the nine points on the panel. 


The procedure is 


to minimize the error function E. 


g 

E - r 2 [ 6 1 + e 2 q + p 3"j + Mj 2 + 6 5 { j"j 

2 2 
+ 3^ - Mj] 


(79) 


In order to guarantee that the panel center doublet strength is 
obtained, 3-i is set exactly equal to yg. The weighting factor 

Wj^ is equal to 1.00, unless j is an edge control point, in which 

case the value 10 + ^ is used to force the doublet distribution 
through the control point value. If point j is an edge control 
point, the values g j and nj are replaced by (1 - e)€j and 
(1 - e)nj, where e (<<1) is an input value that dictates the 
retraction of edge control points. Equation (79) is minimized 
with respect to the five unknowns ^4' ^5' an ^ ^6' This 

leads to a five by five system of linear equations. Matrix in- 

version and multiplication then provide the values b. . for 
equation (78) . 


In principle, it would be possible to place nine boundary 
condition control points on each body panel, where each control 
point would correspond to one unknown value y j . However, the 
computing time would be prohibitive. It is more desirable to 
maintain the control points of figure 26 and establish each yj 
as a weighted average of the doublet strengths at neighboring 
control point locations. For wake sections, the weighted aver- 
age is established under the constraint that the doublet strength 
be constant along the panel edges aligned with any fixed value 
of i in figure 25. By substituting the weighted average into 
equation (78) and then substituting equation (78) into (77), the 
quadratic doublet distribution on a panel is determined as a 
function of the unknown control point doublet strengths. The 
procedure for establishing the weighted average is the subject 
of the following section. 


Doublet Strength at Panel Edges 


Illustrated in figure 28 are the control stations that 
determine the quadratic doublet distributions of all panels in 
the section. In other words, if the doublet strength at each 
point (ip, jy) were known, then equation (78) would provide 
quadratic coefficients for equation (77) . The fact that adjacent 
panels share common control stations at corners and edge midpoints 
limits the magnitude of doublet discontinuity along panel edges, 
regardless of panel spacing. Only the control stations that 
are also boundary condition control points are independent. 
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The doublet strength at each of the remaining stations is deter- 
mined by a weighted average of the strengths at neighboring 
control points. This section describes the approach by which 
that weighting is established. 


(2n — 1 ) 



Figure 28. Doublet Spline Control Stations on a Section 


Cumulative indexing is defined for the doublet control 
stations of all sections. 


k = k (i ) + (i - 1) • (2n-l) + j (80) 

y ^O S y y 

where k,. (i^) is the total number of control stations on pre- 
h*o S 

ceding sections 1 through (ig - 1) . The doublet strength for 

each k y will be a weighted average of the doublet strengths at n D 

neighboring k 's, where each neighboring k must be a control u 

point. The neighboring locations are identified by index ky^, 

where i=l, 2, 3, ..., np. The objective is to determine 

the stations k,. . and weightings c. such that 
** 1 1 


y(k y ) 



E c. * p(k„.) 
i=l 


(81) 


If ky is a control point, the weighting is trivial; i.e., 
n Q = 1, ky^ = k^ , and c^ = 1. In general, quadratic interpola- 
tion is required to determine the proper weighting. The inter- 
polation technique will be described for panels on a body section. 
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Figures 29 and 30 are schematics of neighboring control 
points (1 <_ i <_ n^) surrounding control station k y on a body 
section. Figure 29 applies if ky is a corner point, i.e., if 
the values ip and jp are both odd. Figure 30 applies if ip is 
even and jp is odd. Interchanging ip and jp in figure 30 corres- 
ponds to odd iy and even j . If the separation between ky and 
every edge of the section on which k p lies is no less than two 
panels, then all of the iiq neighboring control points are at panel 
centers. Otherwise, some of the n^ points will be edge control 
points of the section (corresponding to ETYPA = 1 ) or control 
points on the opposite side of the section (corresponding to 
ETYPA ^ 1) . In the former case, the indexing is handled auto- 
matically, and, if the point k y is sufficiently close to the 
section edge, a single edge control point can be identified by 
more than one value i. In the latter case, it is the respon- 
sibility of the user to select the values ky^ corresponding to 

neighboring control points on the opposite side of the section 
edge. This will allow the program to carry the quadratic inter- 
polation across the section edge. If the edge lies on a plane 
of symmetry, the mirror image of point k is assigned 

the same control station index k . p 

u 
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Figure 29. Local Indexing at a Doublet Corner Control Station 


59 


i m mi i i 




II ■ I I II 1 1 II I I II I 


III 




• Dependent Point (i y , jy) 



Figure 30. Local Indexing at a Doublet Edge Midpoint Control Station 


If one or more wake section edges intersect any of the body 
panel edges within the space occupied by the n^ control points, 
the least squares fit will not be applied to the doublet distri- 
bution vi, but to a function u illustrated in figure 31 . In the 
vicinity of point k y , the function u is equal to the doublet 
density p. However, on the opposite side of each wake inter- 
section, u is adjusted by +Ap, where Ap is the discontinuity 
in doublet density on the body surface generated by the wake. 
Therefore, the function u is continuous with continuous gradient 
and is suitable for quadratic interpolation. 

At this stage, it is desirable to establish a local (£,n,C)- 
coordinate system such that the origin is at k y and_^tlie 5-axis 
is locally normal to the body surface. The vector P. is defined 
for 1 _< i <_ n . 


P . = 

1 


where the positive sign for y and z is selected unless k y . is 
on the opposite side of the y = 0 and z = 0 symmetry plants, 
respectively. Vector is the displacement between points k y . 
and k y . The unit vector e is a weighted average of local panel 
normal vectors. ^ 
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Figure 31. Doublet Discontinuity at Wake-Fuselage Intersection 


-> 



c 


n 


D 


{n , + n , + n }, , . . 

x' - y' - z k p (l) 


i ' 1 fTi 


2 2 T 

+ P . 0 + P. _ 

i2 i3 


(83) 


where the constant c accounts for the fact that the right-hand- 

side is not normalized, where kp(i) is the panel on which point 

k,, . lies, where (n , n , n ) is the unit normal vector of panel 
yj_ • x y z 

kp(i) in the general coordinate system, and where the summa- 
tion is to include only values i that are panel centers less 
than one panel distant from k y . Once again, the sign selection 
depends upon whether point i is on the opposite side of a sym- 
metry plane, e^ is defined by selecting any unit vector normal 
to e^ . It follows that e n is defined as e^ x e^- Each point i 
is then assigned coordinates in the (?, n , c ) -system. 
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' 5 i 


= S. . e. 


— fr- 

n . = P. 
x x 


(1 £ i 1 n D ) 


(84) 


= P i ’ e ? 

In order to approximately account for surface curvature, ?p and 
n-j should both be magnified by the factor 
[U 2 + n? + ^t)/(5 2 + t 2 )] 1/2 

X X X X X 

Function u(£, ti ) is modeled as a quadratic 

u(5, n) = 3 ] _ + 3 2 ? + 6 3 n + B + B 5 ?n + 6 6 n 2 (85) 

Designating up as the value of u at point (1 < i < n D ) , 

quadratic interpolation is performed by minimizing function E. 
n D 

E . £ w. 2 (B + qg + q-y + qs 2 + qqq + qq 2 - q> 

X=1 


( 86 ) 


+ 10 4 (B 4 2 + 3 5 2 + 6 6 Z ) 


2 

where Wp is a weighting factor inversely proportional to the 
fourth power of distance between kyp and ky . If ky p is less than 
one panel width from k p , Wp 2 is increased by a factor of 10 +4 to 
force a close fit. The term 


-4 2 2 2 

io (s 4 z + b 5 + b 6 ; 


eliminates the possibility of unbounded curvature in the quadra- 
tic. E is minimized with respect to the unknowns Bp, . . ., $6« 

The resulting system of linear equations is manipulated by 
standard matrix inversion and matrix multiplication techniques 
to provide the array c^ where 


B 


1 


E 

i=l 


c . u . 

X X 


(87) 


The array c^ should be normalized such that E c. is equal to 

i=l 

one. Since the origin of the (£, n, c) -system is at k^ , equation 
(85) implies that 


62 



( 88 ) 


y (k y ) = e-j_ 


Combining equations (87) and (88), 


y (k ) = E c . u . (89) 

y i=l 1 1 

Equation (89) can be transformed into the desired form of 
equation (81) by expressing u^ as y (k p ^) + Ay^ + Ay 2 + . . 

where each Ayj corresponds to a wake intersection somewhere be- 
tween points k y and i. As stated earlier, only intersections 
along wake edge numbers ig E = 2 or 4 are permitted. The doublet 
strength along the wake edge is constant, equal to the doublet 
strength at one of the two corner control points of wake edge 
ig E = 1. Therefore, each Ayj can be replaced by a wake control 
point doublet strength. Equation (89) becomes 


y (k ) 
y 


(n D + V 


E 

i=l 


c . 
i 


y (k ) 


(90) 


where i > n D refers to wake control points. 

If k y is a panel corner on a body or wake section edge of 
ETYPA = 1, the quadratic interpolation is modified to be based 
upon distance along the edge. The modified interpolation in- 
volves a single variable, as opposed to the two variables £ and 
n required for the surface interpolation. For a wake section, 
y (k ) is constrained to be a constant for all values j y at a fixed 
i . This constraint corresponds to the approximation of zero 
pressure loading on a wake. 

By substituting equation (90) into equation (78) , the 
quadratic doublet coefficients on a panel are expressed in terms 
of the independent unknowns at control points. 


B . = E b. . [ E c • y (k ) ] (91) 

1 j = l ^ £ = 1 £ 


Boundary Condition Satisfaction 

There are two types of boundary conditions to be satisfied. 
At each body control point, the internal perturbation potential 
induced by the simultaneous action of all singularities is set 
equal to zero. At each wake control point, the average total 
velocity between the adjacent upper and lower Kutta condition 
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points is to be zero. This average velocity can be expressed as 
the difference in potential between the two points, less the 
local doublet strength on the wake at the control point. 

For each control point there is one boundary condition to be 
satisfied and one unknown doublet strength. The boundary con- 
dition equations are established by using equations (76) and (91) 
in conjunction with available potential influence functions for 
unit distributions of constant source, constant doublet, linear 
doublet, and second order doublet density (reference 18). 

The resulting system of linear equations is solved for the 
doublet strength at each control point. Equations (77) and (91) 
then provide the doublet distribution on every panel. 

Velocity and Pressure Calculation 

The calculated doublet distribution is converted to surface 
velocity and pressure at each body panel center, and then the 
forces and moments are integrated. For each body panel, the 
flow velocity in the local panel coordinate system is 


V 

= v -e. 

+ 

8 

€ 

“ 5 




->■ 



V 

= v -e 

+ 

8 

n 

co r| 




->- -b 



V 

= V *e 

+ 

a 

e 

CO £ 




(92) 


Equation (92) expresses the fact that on a body surface the local 
source strength is equal to the perturbation velocity normal 
component and the local doublet strength is the local perturba- 
tion potential, corresponding to Green's third identity. For 
steady flow, Bernoulli’s equation provides the local pressure 
coefficient . 


C 

P 


2 2 2 
v/ + V + V 

K n Z 


V 


ref 


(93) 


Forces and moments are calculated under the approximations 
that the pressure is uniform on each panel and that the flat 
quadrilateral panel is the local geometry. 

This completes the description of the analysis method of 
the MCAIR 3-D Program. 
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Example Solution 


The MCAIR 3-D Program has been checked out for a variety 
of simple geometries. At the time of writing this report, the 
program had not been evaluated for wing-fuselage configurations. 
Although it is anticipated that the evaluation will be success- 
ful, only a simple isolated wing solution is available for pre- 
sentation here. 

The right half-wing of figure 32 was represented by five 
spanwise strips of panels with 10 upper and 10 lower surface 
panels per strip. Twenty additional panels were used to close 
the tip. The wake was modelled as being tangent to the local 
trailing edge bisector angle at each span station. 


The calculated incompressible inviscid pressure distribu- 
tion at two angles of attack is shown in figures 33 and 34. 
Unpublished experimental data is also shown for comparison. The 
data was taken in the presence of a fuselage at 0.6 Mach number 
and 6 million Reynolds number, based on the mean aerodynamic 
chord. Although no corrections were made for fuselage inter- 
ference, compressibility, or viscosity, the calculated results 
agree reasonably well with the experimental data. The solution 
required approximately 50 seconds computing time on the CDC CYBER 


175. 


FIRST-ORDER EXPANSION APPROACH 
FOR THICK WING-FUSELAGES 

The most widely recognized feature of three-dimensional 
surface singularity panel methods is that the numerical solution 
formulation is not dependent upon the type of surface geometry. 

A second feature that has been almost totally unrecognized is 
that the analytical expressions in a panel method formulation 
can be linearized to provide the matrix of velocity derivatives 
with respect to surface geometry perturbations. This approach 
was used for the design mode of the Multi-Element Airfoil 
Inviscid Analysis and Design Program (MAAD) , and it is antici- 
pated that the success of the two-dimensional design method will 
carry over to three-dimensional configurations. 

Possibly an even more valuable application of the velocity- 
geometry derivative matrix is to the analysis of several wing- 
fuselage geometry variations at no significant additional ex- 
pense. During the configuration development of most airplanes, 
several possible combinations of wing planform, twist, and 
section shape are considered. With existing surface panel methods, 
each combination would require one independent analysis solution. 

If compressibility effects were to be included, each Mach number 
would also require a separate solution. Since a three-dimensional 
analysis can cost several hundred dollars and require overnight 
turnaround, there is reluctance among many engineers to exercise 
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Figure 33. WLB-1 Pressure Distribution at 3.8° a 













surface panel methods. On the other hand, 
plane panel methods such as vortex lattice 
ly accurate, particularly at leading edges 
intersections. It is anticipated that the 
viated by implementing a velocity-geometry 
lation in the MCAIR 3-D program. 


the simpler chord 
are often insuf f icient- 
and wing-fuselage 
dilemma would be alle- 
perturbation formu- 


Outlined in this section is the approach through which 
arbitrary small variations in both configuration geometry and 
subsonic Mach number can be analyzed at very small cost. Then 
an approach for solving finite wing design problems will be 
summarized. 


Low Cost Method for Analyzing 
Geometry Perturbations 

Consider any three-dimensional geometry subject to analysis 
by the MCAIR 3-D program. Suppose that during the analysis, 
a velocity-geometry derivative matrix were calculated and stored 
on a computer file. At any future date, it would be possible 
to predict the effect of a perturbation in configuration geometry 
or subsonic Mach number by retrieving the matrix and multiplying 
it by a new right-hand-side. The additional cost would be insig- 
nificant. Applications would include configuration evaluation 
and optimization, viscous-inviscid interactions, and the calcu- 
lation of stability derivatives. 


The discussion to follow will concentrate on the utility of 
a precalculated velocity-geometry derivative matrix, not on 
the method used to establish the matrix. That method has not 
yet been formulated for the MCAIR 3-D program. However, it is 
expected to be a straightforward extension of the existing 2-D 
method incorporated in Program MAAD . In order to simplify the 
discussion, only zero normal velocity boundary conditions will 
be considered. It is understood that a general perturbation 
solution formulation can and indeed should include the option 
for nonzero normal velocity. Attention will initially focus on 
simple two-dimensional flow past an airfoil. Then three- 
dimensional geometry perturbations and Mach number corrections 
will be discussed. 

Consider the airfoil of figure 35. The baseline geometry 
is perturbed by displacing each panel endpoint j by . Al- 
though the panel indexing remains fixed, the coordinates are 
translated. Arbitrary perturbations are permitted in airfoil 
chord, panel length distribution, and airfoil location with res- 
pect to any other elements in the flow. The velocity at the 
midpoint of panel i before and after the geometry perturbation 
is respectively designated Vj_ and (Vj_ + 6Vp) . Assume that the 
velocity V<x> is parallel to the x-axis and of unit 
i.e. V m = e x . Then a velocity-geometry derivative 
can be established such that 


f reestream 
magnitude, 
matrix X . . 
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V: + 5 V; 

* 1 Geometry 



Figure 35. Geometry - Velocity Perturbation for an Airfoil 


6V = Z (A -SR ) + 0(|6R.[ 2 ) (94) 

1 -*-l J 3 


The vector nomenclature for matrix Aj_j is intended to indicate 
that there are two matrix coefficients for each (i, j) , respec- 
tively corresponding to the x and y components of 6Rj . The 
derivative matrix is calculated as a function orthe baseline 

geometry only. The effect of the geometry perturbation is pre- 
dicted by multiplying Apj and 6Rj . This neglects second order 
and higher terms of 6$-; , identified by the nomenclature 
OllSR-jl 2 ). 3 

Although only linear terms are included, the prediction 
accuracy can be good for surprisingly large perturbations. An 
example is the perturbation from a circular cylinder to a foot- 
ball shape (figure 36). Using Program MAAD , the change in 
velocity distribution was calculated by three approaches. First, 
the actual change was calculated by independently solving for 
the flow about the perturbed and unperturbed ^geometries . Second, 
the present approach of multiplying Apj by 6Rj was applied. 

Third, the method of equivalent blowing, also known as the sur- 
face transpiration method, was used. The present first-order 
method provides substantially better agreement with the actual 
velocity change than the method of equivalent blowing (figure 
36) . The reason is that only in the present approach are all 
first-order perturbations included. 
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Now the procedure for handling arbitrary angles of attack 
will be described. It is necessary to introduce a second 
velocity-geometry derivative matrix, corresponding to V m = ey. 
Superscripts (1) and (2) are respectively used to denote baseline 
flow solutions for V m = e x and = ey. Suppose that the free- 
stream velocity V„> is at an angle of attack a with respect to the 
x-axis. The baseline solution is 

V. = V [V. (1) cosa + V. (2) sina] (95) 

l » l l 


The first-order velocity change induced by the geometry pertur- 
bation at the fixed a is 

<5V. = + V cosa E (i. (1) -<5R.) + V sina E (£.^*s£.) (96) 

1 oo j 11 1 “ j 11 3 


By using equations (95) and (96) , it is possible to isolate the 
effects of geometry perturbation and changes in angle of attack. 


For three-dimensional conf igurations , there are three 
components to the geometry perturbation vector and two com- 

ponents of tangential velocity at each panel center. Consider 
the velocity in the local (e, , rw £ ) -coordinate system of any 
panel. Before the geometry is perturbed, the components are 


(V; 


V 


eral ( x > Y > z ) 


0) . If the free stream velocity components in the gen- 

) , then 


system are (Vco-^/ 

Voo 2 

' Vco 3 

(1) 
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v< 2 > + 
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v 5 (3 > 

(1) 
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V (2 > + 


v (3) 


(97) 


where superscripts (1), (2),_^(3) respectively correspond to 

baseline flow solutions for = e x , = ey, = e z . Following 
the geometry perturbation, the velocity components in the updated 
(?, n, C ) -system will be (V^ + 6V^ , V n + 6V n , 0). The three- 
dimensional analogue to equation (96) can be expressed in the 
form 
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( 98 ) 


72 


where the j -summation will include only a few neighboring corner 
points and the k-summation will include every corner point of 

the panelled configuration. The coefficients , a^), h r . , b . 

Sd in Zj' no 

and bK ' are functions of only the baseline geometry. The signi- 

^ , -v ( 1 _ \ 

fxcant computing expense is associated with generating Bj k . 
Usually, only values for i = 1 and 3 will be calculated, with 
the x-z plane being a plane of flow symmetry.^ For a configura- 
tion modelled by one thousand panels, matrix B^) will contain 
approximately three million values for each of J the three values 
of i. It is anticipated that the calculation of the coefficients 
of equation (98) will increase the expense of a conventional 
analysis solution by a factor of five. The payoff is that once 
the coefficients are generated, equation (98) will provide 
rapid, inexpensive predictions for arbitrary geometry perturba- 
tions. Approximate force and moment perturbations can be calcu- 
lated by incorporating equation (98) in the Bernoulli equation 
and then integrating the resulting pressures over the perturbed 
geometry. 


Mach number variations can be simulated by using equation 
(98) to predict the effect of the geometry stretching associated 
with Gothert's rule. For example, suppose that at a given 
baseline subsonic Mach number and angle of attack, the equivalent 
incompressible coordinates of some body point are (x, y, z) . 

At some other free stream conditions, the equivalent incompressi- 
ble coordinates are (x, y, z) , where 


x | 

! 

b ll b 12 b 13 

lx 

1 >1 


b 21 b 22 b 23 

Jy 

1 N 


b 31 b 32 b 33_ 

1 z 


The 3x3 matrix b^j is established by Gothert's rule. Equation 
(99) can be written in the form 


x = x + <5x 

y = y + <5y (100) 

z = z + (5 z 


where 6x, <5y , and Sz are functions of (x, y, z) . Equation (100) 
is then substituted directly into equation (98) to predict the 
effect of a change in freestream Mach number. 

The benefits of the present perturbation analysis approach 
will have the greatest impact when used in conjunction with 
suitable geometry and interactive graphics packages. For example. 
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suppose the user digitizes a new set of planform coordinates; 
the geometry package would generate the corresponding values 
for <5Rj, retrieve the velocity-geometry derivative matrix, and 
perform the required calculations; and the interactive graphics 
facility would display the desired pressures, forces, and 
moments. The entire process would require only a few minutes at 
no appreciable expense. 

Finite Wing Design 

In an earlier section, an iterative method was presented 
for calculating the geometry corresponding to a prescribed air- 
foil pressure distribution. In each iteration cycle, the geome- 
try would be updated and a new analysis solution calculated. 

This procedure is reasonable for two-dimensional flow because 
the associated computing cost is usually quite small. However, 
for three-dimensional configurations, an analogous design method 
could be prohibitively expensive. Described here is an approach 
that should provide a more economical method for designing the 
section geometries of a finite wing. 

The approach involves iteration, using a given starting 
geometry to initialize the calculations. Unlike the two- 
dimensional design method, the starting geometry will not be up- 
dated during every iteration cycle. Significant computing time 
will be saved because the baseline analysis solution and velocity- 
geometry derivative matrix can be reused for several iterations. 
Then the starting geometry can be updated and the process re- 
peated. A flow chart of the design approach is presented in 
figure 37. Each step is described below. 

The objective is to determine the wing streamwise section 
geometries at fixed spanwise stations corresponding to a pre- 
scribed pressure or velocity distribution. The fuselage geometry 
and wing planform are assumed to be known beforehand and will 
not change. If a velocity distribution is prescribed, both the 
streamwise and spanwise components must be specified. If a 
pressure distribution is prescribed, a wing internal volume con- 
straint should be applied to eliminate the possibility of a 
nonunique solution. To enhance numerical stability, prescribed 
values will apply to the centers and edges of the wing 
panels . 

The MCAIR 3-D Program will be used to analyze the starting 
geometry. The analysis will establish the difference between 
the initial and the desired distributions of pressure (or velo- 
city) . Then the velocity-geometry derivative matrix of equation 
(98) will be calculated. 

Equation (98) will be inverted to calculate the geometry 
perturbation array 6Rj associated with the prescribed pressure 
distribution. The inversion will be established in the sense 
of the least-square-error at panel centers, corners, and edge mid- 
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Figure 37. Design Approach for Finite Wings 

points. Appropriate linear constraints should be applied for 
geometric closure and, possibly, prescribed internal volume. 

Also, the length of the panel streamwise edges should be frozen 
during the inversion in order to prevent unnecessary disturbance 
to the surface distance distribution of the panelling. If the 
calculated array <5Rj is unrealistic or insufficiently smooth, 
the values can thenbe modified manually. At that stage, it will 
not be necessary to maintain the same lengths of panel streamwise 
edges as in the starting geometry. 

By applying equation (98) and Bernoulli's equation, the cor- 
rected pressure distribution is calculated, corresponding to 
the starting geometry perturbed by . It is noted that the 
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corrected pressure distribution (C p + 6C p ) is only accurate to 
first-order with respect to SRj . If the p dif f erence between the 
corrected and prescribed pressure distributions is unusually 
large, it is possible that the prescribed distribution does not 
exist. In that case, the corrected pressure distribution should 
guide the user in determining the smallest modification required 
to create an existing prescribed distribution. If the agreement 
between the prescribed and corrected pressure distributions is 
sufficiently good, the array 6Rj will be used to update the 
starting geometry and thg entire procedure will be repeated. 
Otherwise, a new array <5Rj should be calculated by again invert- 
ing equation (98). 

The present finite wing design approach incorporates several 
features that proved essential to two-dimensional design. Based 
upon experience with two-dimensional design, it is expected 
that no more than one or two starting geometry updates will be 
required to obtain a converged three-dimensional design solution. 


CONCLUSIONS 

The use of combined source-doublet panels leads to accurate 
and numerically stable subsonic analysis and design solutions. 

The prediction accuracy associated with internal perturbation 
potential boundary conditions and flat panels is competitive 
with more complex curved panel formulations. 

The Multi-Element Airfoil Inviscid Analysis and Design 
Program can calculate the surface shape and arclength in design 
regions, corresponding to a prescribed pressure distribution. 
Converged solutions are usually obtained in fewer than five 
iteration cycles. The rapid convergence results from the use of 
all first order velocity-geometry perturbation terms in the in- 
verse step of each iteration cycle. The variable arclength 
feature can automatically transform a nonexisting prescribed 
pressure distribution into one of the closest possible distribu- 
tions, thereby simplifying user input requirements. Applications 
of the design option include airfoil optimization and viscous- 
inviscid interactions. 

The MCAIR 3-D Subsonic Potential Flow Program predicts the 
inviscid pressure distribution for arbitrary lifting or non- 
lifting configurations. The method of least squares minimizes 
doublet discontinuities at panel edges, which in turn reduces 
the effect of the selected panelling arrangement on the numeri- 
cal solution. It is possible to perform a first-order expansion 
to the solution formulation of the program. This would allow 
the calculation of a velocity-geometry derivative matrix that 
could be stored for future use. The analysis of an arbitrary 
configuration modification would simply require an inexpensive 
multiplication of the matrix by a new right-hand-side. The 
matrix would also provide the foundation for a three-dimensional 
wing design method. 
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